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Abstract 

A new method is proposed for integrating the equations of motion of an elastic filament. In the 
standard finite-difference and finite-element formulations the continuum equations of motion are 
discretized in space and time, but it is then difficult to ensure that the Hamiltonian structure of the 
exact equations is preserved. Here we discretize the Hamiltonian itself, expressed as a line integral 
over the contour of the filament. This discrete representation of the continuum filament can then 
be integrated by one of the explicit symplectic integrators frequently used in molecular dynamics. 
The model systematically approximates the continuum partial differential equations, but has the 
same level of computational complexity as molecular dynamics and is constraint free. Numerical 
tests show that the algorithm is much more stable than a finite-difference formulation and can be 
used for high aspect ratio filaments, such as actin. 
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I. INTRODUCTION 



Elastic rods are a ubiquitous model of semi-flexible biopolymers such as DNA,-'^'^'^'^'^'^ 
actin,-ii^iiiii^ and microtubules.— They can also be found in a diverse range of applica- 
tions including catheter navigation,— undersea cables,— and organismal biology.— In bio- 
physics, the worm-like chain (WLC) model^*^ underpins many theoretical^»i^'i^'2Si^iS2ii^;^ 
and numerical^'2^'2Ii2^ studies of semiflexible polymers. The WLC model is a linearization 
of the classical Kirchoff rod model,^'^ which is itself a limiting case where the product of 
the local curvature and filament thickness is everywhere small.— In this limit the shear and 
extensional strains are negligible but the constraint forces generated by them are not. In 
this paper we consider a generalization of the Kirchoff model,— where the shear and exten- 
sional strains are explicitly accounted for by an elastic constitutive model, eliminating the 
need for constraint forces at the cost of an additional time scale; such models are frequently 
referred to as "geometrically exact" in the finite-element literature.— 

The dynamics of Kirchoff or geometrically exact (GE) filaments is typically determined by 
finite-element or finite-difference approximations, but the stiffness of the numerical system 
has proved to be a difficult and long-standing problem.— Significant progress has been 
made by developing implicit methods that exactly satisfy the constraints of momentum and 
energy conservation,—!^ yet even here artificial dissipation is often needed for long-term 
stability.^'' On the other hand, in discrete dynamical systems it is known that symplectic 
integration methods give superior long-term stability in comparison with either high-order 
explicit or implicit integration methods;— the most common symplectic integrator is the 
Verlet algorithm.— Symplectic integrators generate a sequence of canonical transformations, 
which do not exactly conserve energy but do preserve the density of points in the phase 
space, along with the Poincare invariants. In recent years symplectic integrators have been 
developed for both linear and angular motions.— i^i^ The objective of this paper is to explore 
a symplectic integration method for geometrically exact filament models. This requires both 
a Hamiltonian approximation to the partial differential equations describing the filament 
dynamics, and a symplectic integrator. 

The proposed algorithm is based on a discretization of the Hamiltonian line integral of an 
elastic filament, including shear and extensional degrees of freedom. Since the nodal forces 
and torques follow from an exact differentiation of a potential function, the equations of 
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motion are guaranteed to be Hamiltonian, although the potential function itself is only an 
approximation to the continuum limit. This is in contrast to finite-element methods, where 
the continuum equations of motion are discretized in space; in this case the Hamiltonian 
structure is not preserved, even if the total energy is conserved.— In fact, it can be shown 
that for any approximate solution it is not possible to maintain both the symplectic structure 
and exact energy conservation simultaneously.— 

An outline of the paper is as follows. In Sec. [TTl we describe different models of elastic 
filaments-GE, Kirchoff, WLC-and indicate how they are related. Next (Sec. IIII|) . we derive 
a simple finite-difference approximation of the equations of motion of a GE filament model, 
as a basis for comparison with the Hamiltonian formulation presented in Sec. IIVI We note 
that the Hamiltonian approach has only been followed occasionally,^^ and in that case for the 
Kirchoff rod model. We will argue (Sec. |V]) that the absence of geometric constraints in the 
GE model offers computational advantages over the Kirchoff model when there are excluded 
volume interactions between the segments. We replace the usual implicit time integration^^i^ 
with an explicit operator splitting method,"^ which eliminates the repeated force evaluations 
of an implicit method. The numerical scheme is stable and energy conserving even for large 
deformations; we illustrate this by numerical example in Sec. |Vl Our conclusions and future 
outlook are in Sec. IVII 

II. ELASTIC FILAMENT MODELS 

The classical Kirchoff theory of elastic rods has been elegantly and concisely described in 
the " Theory of Elasticiti/^ by Landau and Lifshitz,— and the seminal book by Love.— More 
rigorous derivations of the equations of motion are available in the literature.— >^ Here we 
summarize the key concepts and establish the notation to be used later in the paper. An 
elastic filament (or thin rod) is described by the coordinates of its centerline r{s) and a set 
of orthonormal directors di{s), ^2(5), 0^3(5). The directors establish the orientation of a 
cross section or material plane at the location s, where s is a parametric coordinate defining 
the position of each point along the centerline. In the undeformed filament, s is the contour 
length from the origin. We will choose a body-fixed coordinate system such that di and d2 
point along the principal axes of inertia of the cross section and therefore ^3 = di x ^2 is 
normal to the material plane; the coordinate system is illustrated in Fig [H If the rod has 
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FIG. 1: An elastic filament in the unstrained (reference) state (a) and after deformation (b). In the 
reference state, the material plane, shown by the solid ellipse, is aligned with its normal parallel to 
the tangent to the centerline (dashed line). The local director basis of the reference state, d^{s), 
and the deformed state, di{s), are also shown. A material point (solid black circle) moves with the 
translation and rotation of the local coordinate system; in this case extension, shear, bend, and 
twist can all be seen. 

a circular cross section then the initial choice of di and d2 contains an arbitrary rotation 
about 0^3. In contrast with the Kirchoff theory, we will not assume that ^3 is constrained to 
be parallel to the tangent vector dgr (Fig[T)D). 

The key assumption of thin-rod elasticity is that there is no deformation within a material 
plane, only translation and rotation of that plane. Deformation of an elastic filament is then 
described by two one-dimensional strain fields, r(s) and 0(s), describing the rate of change 



4 



of the centerline position and director vectors along the filament^i^^ 

Ti = di ■ (dsv) ^1 = ^3- {dsd2) = -d2 ■ (S^dg) 

T2 = d2- (dsv) ^2 = di ■ (dsd^) = -d3 ■ (dsdi) (1) 

Ts = ds- (dsr) ^3 = da ■ (dsdi) = -di ■ (S.da). 

A thin segment of the filament can be subjected to six different deformations. Fi and r2 
describe transverse motions of a material plane with respect to the normal vector (^3), 
which causes shearing of the segment, while Fs describes extension or compression of the 
segment. Bending of the segment about its principal axes is described by Vti and ^2, and 
twisting of the segment by ^3. Uniform deformation corresponds to constant values of T 
and O; for example, in a straight rod T = [0,0,1] and Q, = [0,0,0]. More interestingly, a 
helical rod can be described by a constant bend and twist, T = [0,0, 1], ft = [Rn^ ,0, Pk^], 
where R is the radius of the helix, 2ttP is the pitch, and the combined curvature due 
to bend and twist, k = (P^ + The choice of signs define a right-handed helix, 

r{s) = [-Rcos(ks), -Rsin(/ts), P/ts], with basis vectors 

di = [P/t sin(Ks), — Pk cos(ks), Pk] 

^2 = [cos(ks), sin(Ks), 0] (2) 
^3 = [— PKsin(Ks), Pkcos(ks), Pk]. 

The stresses in the rod are assumed to be linear in the deviations in the strain fields, 
AFj = Fj — F? and Af2j = fij — from the reference (stress free) configuration r°, 0°. 
It is convenient to define the strains in the body-fixed coordinate system, since the elastic 
constant matrix is then diagonal. The force F[ and couple Fp on each material plane are^^i^ 

F[ = ClAT,, FP = CfAn,, (3) 

where the elastic constants for each deformation are, in principle, independent. In the GE 
model, the strain energy density U{s) contains contributions from shear and extension, in 
addition to the usual bend and twist of the Kirchoff model, 

u = u^ + u^ = \Y. {cl^n + cpAn^) . (4) 

i=l 

For an isotropic material, the elastic moduli for shear (Cfg), extension (C3 ), bend (CP2)) 
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and twist (C^) are given by: 

Cl = GA = Yh 

Cl = CA C^ = Yh (5) 

Cl = YA C^ = Gh, 

where G is the shear modulus, Y is Young's modulus, A is the area of the cross-section 
and Ji and I2 are its principle moments of inertia. For rods with a circular cross section, 
J3 = Ji + /2, but in the general case there is an additional contribution from the warping of 
the cross section, so that J3 is then distinct from Ii + h- The elastic coefficients can also 
be determined empirically, without reference to any particular constitutive law. 

The velocity and angular velocity of the segment are defined in an analogous fashion to 
the strain fields in Eq. [H 

vi = di- (dtr) uji = d3- {dtdi) = -^2 • (dtds) 

V2 = d2- (dtr) uj2 = di- {dtd^) = -dg • (dtdi) (6) 

^3 = • (dtr) = da ■ (dtdi) = -di ■ (Sfds). 

The kinetic energy density of the filament is then^^i^ 

1 ^ 

T = + = - J] {M[v^ + MfuJ) , (7) 

i=l 

where the generalized mass densities associated with shear (Mf , Ml"), extension (M3 ), bend 
(Mf ,M2^) and twist (Mg^), are 

Mf = pA, MP = p/„ (8) 

and p is the mass density of the filament. 

Equations of motion for the filament can be derived from the balance of linear and angular 
momenta in a thin segment bounded by the planes s and s + ds. The rate of change of the 
linear momentum of the segment, pds, is 

pds = F^{s + ds) - F^{s), (9) 

where p = ^^=1 M'^fjdj is the linear momentum density (per unit length). The forces on 
the two planes must be differenced in a common coordinate frame, which we take as the 
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space-fixed frame. Tlie balance of angular momentum in the segment Ids involves both 
couples and moments of the force, 

ids = F^{s + ds) - F^{s) + r{s + ds) x F^{s + ds) - r{s) x F^{s), (10) 

where I = X]^=i MpUidi is the linear angular momentum density. Thus the equations of 
motion of a GE filament are 

P = dsF^, (11) 

i = dsF^ + r' X F^, (12) 

where r' = dsV indicates a spatial derivative along the filament. A finite-difference approx- 
imation to these equations is described in Sec. Illli 

Equations [TTHT2] describe the dynamics of the GE rod model.— The difference with the 
Kirchoff theory is that, here, the force on a material plane, F[, is given by a constitutive 
equation, Eq. ([3]), based on the deflection and extension of the local tangent vector relative 
to the material plane, Eq ([T]). In the Kirchoff model the tangent vector is constrained to 
remain parallel to (unshearable) and of unit length (inextensible), or in other words 
AFj = and r' = d^. As a result, neighboring segments can only rotate with respect to one 
another, leading to a compatibility condition,— 

v' = ujxr' = dg, (13) 

where the last equality follows from the kinematic conditions, di = uj x di.—^ Differen- 
tiating Eq. ffTTl) with respect to s gives an equation for the constraint force satisfying the 
compatibihty equation, 

d^,F^ = M^ds, (14) 

where d^ = u: x d^ + uj x [uj x ^3).— i^i^ A simpler, but approximate solution is to neglect 
the angular momentum perpendicular to the tangent vector,—!^ and determine the shear 
forces, i^^'"*", directly from the cross product of Eq. ( JT2l) with 0^3, 

d3 X dsF"" = (1 - ^3^3) ■ = F^'^. (15) 

The force along d^ is determined from the inextensibility condition,— 

dgr ■ dgV = 1. (16) 
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The Kirchoff model has the computational advantage that the shear and extensional modes 
are frozen by the constraints, so that a larger time step may be used. On the other hand 
the numerical integration is inherently implicit and must be solved iteratively at each time 
step. 

Bending forces can also be determined from the curvature in the centerline position 
vector, r' x r", rather than from derivatives of the basis vectors, Eq. ([1]). In the case of 
a weakly bent filament, the tangent can be assumed to be locally constant,— and, with an 
isotropic bending stiffness = = C^, 

Differentiating once more (again ignoring derivatives of r') , we obtain the equation of motion 
for the bending of a WLC,^'^'^^ 

M^r = -C^ (1 - rV) ■ r"", (18) 

although what is really being calculated is the constraint force needed to resist the shear 
deformations arising from the compatibility condition, Eq. [131 In addition, a constraint force 
is needed to satisfy the inextensibility condition, Eq. (I16p . Unfortunately, Eq. [TS] is very 
stiff, and numerical integration of the partial differential equations is not straightforward.— 
Most simulations of the WLC model have therefore discretized the filament into a sequence 
of beads interacting via a bending potential.— i^^i^ Although this sacrifices fidelity to the 
continuum filament model, the ordinary differential equations for the bead positions can 
be integrated using standard molecular dynamics methods, including constraint forces to 
maintain a discrete approximation to Eq. fll6p . In this paper we derive a discrete Hamil- 
tonian representation of a GE rod model, along the lines already established for the WLC. 
Our algorithm systematically approximates the GE filament model, while maintaining the 
simplicity of the WLC approach. We wish to emphasize that the models described in this 
work are discrete approximations to continuous filaments, in which the nodes indicate rep- 
resentative points along the centerline. This is different from models where the segments 
are physical objects with finite length, undergoing rigid-body motion.— 
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III. DISCRETE EQUATIONS OF MOTION 



We first describe a spatial discretization of tlie equations of motion of a GE rod, Eqs. \TT\ - 
rC2[ The filament is divided into equal segments of length As = L/N, and nodes are 
defined at the center of each segment,— 

Sn= {n-D As, n = l,2,...N. (19) 

The instantaneous state of the filament is then given by the nodal coordinates r", quaternions 
g^, linear momenta p^, and angular momenta We use Greek subscripts, 7, to 
indicate components in the space-fixed frame, subscripts i,j,k, to indicate components in 
the body-fixed frame, and the subscripts a, b, c, to denote the components of the quaternion, 
la = [(lo,<lx,(ly,<lz]- The Einstein summation convention is applied to the subscripts 7 
and a, b, c, but not to the indexes i,j, k. Thus for example 

3 

Pa = ^Pidia, Pi = diaPa- (20) 

1=1 

The quaternion Z = [go, q] describes a rotation about an axis parallel to the vector q = 
[Qx, %, Iz] by an angle {} = 2 cos~^(go)- The orientation of a body in space can be specified by 
the components of Z, which we denote by g^. We use quaternions in preference to the director 
basis vectors as angular coordinates,-^ since it reduces the number of degrees of freedom. 
Symplectic integration algorithms using operator splitting exist for both quaternions^'' and 
director vectors. The choice of the body-fixed angular momenta is guided by the integration 
algorithm,— which requires them for the quaternion update. Key properties of quaternions 
are summarized in Table [T] and derived in Appendix |Xl 

An infinitesimal rotation about the body-fixed axes can be written in terms of variations 
in the quaternions (see Appendix |X] for details), 

S(f)i = 2eiaSqa, (21) 

where the quaternion variation is subject to the normalization constraint SqaQa = 0. In 
other words the variation in g^ must be in a three-dimensional space orthogonal to qa- The 
quaternion basis vectors Cj (Eq. IT1.3P describe rotations about a body-fixed axis and are 
orthogonal to each other and to the quaternion itself. The factor of 2 arises because it takes 
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TABLE I: Properties of Quaternions (Appendix lAj) 



qo = cos (I) cos (^^) 
Qr, = sin (I) cos (^^) 
Qy = sin (|) sin (^^) 
= cos (f ) sin (^^) 

Relation between quaternions and Euler angles (0, i/:)^^ 



i52 



(Tl.l) 



^ Qo + qI 



ly - Tz '^{QxQy + qoQz) '^iQxQz - Qoqy) 



\ 



2{qyq.-qoqz) q'o - ql + q'y - q'z 2(g,g, + go?.) 
y 2{q^q^ + qoqy) 2{q^qy - q^q^) ql - ql - ql + ql j 

Director basis in terms of quaternions 



(T1.2) 



62 



\ 



qz -qy 



qo 



\ 



Body-fixed rotations in a quaternion basis 



(T1.3) 



dqa 



— ^ ^ 2tijkdjoL(ika ~\~ 2q(idia. 
j,k=l 

Derivatives of d vectors 



(T1.4) 



dqb 



~ ^ ^ ^ijk^^ja^kb ~l~ ^iaqb qa^ib- 
j,k=l 

Derivatives of e vectors 



(T1.5) 



a product of two quaternions to describe a rotation (Appendix |A]) . The inverse relation 



1 ^ 



(22) 



1=1 



automatically maintains the normalization of g^. The angular velocity and bending strains 
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can be directly related to derivatives of qa, 

Ui = (j)i = 2eiaqa fii = 0- = ^dqq'^. (23) 

We are now in a position to write down ordinary differential equations that approximate 
the dynamics of an elastic filament. A nice feature of the midpoint discretization^'^ is that 
the strains are naturally evaluated at integer multiples of the segment length, nAs, with 
n = 0,1,..., A^. An additional differencing of the internal forces and couples then gives 
accelerations back at the nodal positions. Thus the algorithm is second-order accurate in 
As, with only three nodes directly interacting with one another, just as in the WLC model. 
The derivatives r'^, q'^ are approximated by centered differences at the discrete locations 
nAs, midway between the nodes, 

C = + 0(As)^ (24) 

q': = ^"'''^^^"'^"^^^ + OiAsf. (25) 

In addition we need to estimate the quaternions at nAs in order to calculate the rotation 
matrices, Eqs. IT1.2f[TT73| 

^ = €''\t) + q:it) 

" \qS^'it) + qm\ ^ ^ 

Thus the coordinates, r^,qa, are evaluated at the nodal positions, {n + l/2)As, while the 
derivatives r^", g^'", and mean, q^, are evaluated at riAs. 

The elastic forces and couples at the interior positions nAs, n = 1, 2, . . . , A^ — 1, are then 

3 

Fl'"- = J2 (4^" - r°) , (27) 



i=l 
3 

i=l 



Cfd-, (2e7,gr - ^1) , (28) 



where the notation rf"^ and e"^ indicates the basis vectors are calculated from the average 
quaternions (Eq. [26]) . The forces at the ends of the rod, n = and n = N, are determined 
by the boundary conditions. For free ends, 

F,r,o ^ p^r,N ^ ^n,o ^ pn,N ^ (29) 

while prescribed external forces and couples on the ends of the rod can also be included. 
Dirichlet boundary conditions require virtual nodes, n = and n = A^ + 1, which are 
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constructed to satisfy the boundary conditions at the ends of the filament.— For example, if 
the position and orientation of the rod at s = are specified by f° and g^, then the virtual 
coordinates are 

rl = 2fl - tI (30) 

v/(2g°-gi)(2g2-g^)' 

The elastic forces and couples at s = can then be determined in the same way as for the 
interior nodes. However, it seems preferable to implement Dirichlet conditions by placing 
the nodes at integer locations along the filament, riAs, and then calculating the forces at 
the half-integer positions; this eliminates the need for virtual nodes. In the case of mixed 
boundary conditions a combination of these strategies may be necessary, depending on the 
specifics of the problem; in this paper we just consider filaments with force and couple free 
boundaries. 

The nodal coordinates and momenta satisfy the ordinary differential equations (n = 
1,2,. ..,iV) 

(32) 



pr,n rpV.n—l 

PI = f: = " , (34) 

- *a - I ^ 2^ (^ijkdia ^ I • l-^^j 

\ i,j,k=l J 

The rotation matrices d"^ and e^^, without the overbar (c./. Eqs. |27] and [28]) , are evaluated 
from the nodal quaternions g", whereas the strains F", fi" and forces F['"', F^'"' are evaluated 
at the points nAs, midway between nodes n and n — 1. The numerical approximation to 
the term T x requires nodal values of T and F^ , which are determined by averaging the 
body-fixed strains and forces, and then rotating the vector product to the space-fixed frame 
(Eq.ESD. 



IV. HAMILTONIAN FORMULATION 

The standard procedure for solving the partial differential equations for the linear and 
angular moment a^^'^'^'^'^ does not, in general, lead to a symplectic algorithm, because the 



12 



discrete nodal forces are not derived from a potential energy function. Rather than discretize 
the equations of motion for the continuum rod, we instead discretize the line integral making 
up the Hamiltonian function,'^ to obtain a discrete Hamiltonian that is a second order (in 
As) approximation toTi = T + U. We then use time integration schemes that preserve the 
symplectic structure of the discrete Hamiltonian.— 



A. Hamiltonian for an elastic filament 



The kinetic (Eq. [7j) and potential (Eq. H]) energies of an elastic filament can be written 
in terms of the coordinates and their space and time derivatives, 



r 

2 



^M^f„f„ + 4 J^MfeiaCifegagftj (36) 
^ = \til [Clid.-^ya - ^'){d^pT'^ - r°) + Cf{2e,aq'a - ^^°)(2e.,g^ - fi°)] ds. (37) 

i=l 

The first step is to identify the momentum fields, P = dT/dQ, conjugate to our chosen 
coordinates, Q{s,t) = [ra{s,t),qa{s,t)]: 

3 

Pa = M^r^, la = 4:J2 Mpeiaeibqb, (38) 

i=l 

where la = [lojxjyjz] is the angular momentum field conjugate to qa- It is related to the 
body-fixed angular momentum field, k = MpUi = 2Mpeibqb, 

^- = '^Y.^i^^-^ /. = -e,X. (39) 

i=l 

Rewriting the kinetic energy in terms of the conjugate momenta, 

we can derive the equations of motion of the coordinates by functional differentiation of 
T(P, Q) with respect to P: 

Pa 



3 7 1^7 



i=l ' i=l 
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The equation of motion for the hnear momentum field derives from the potential energy 
due to shear and extension (Eq. [37|) . 



SU^ [•■ ^ir'p 



Flj±ds. (43) 

a Jo 'Ji a 

The functional derivative requires an integration by parts to convert variations in r' to 
variations in r, 

Pa = d,F^, (44) 

as before (Eq. [TTj) . Here we have omitted contributions derived from work done on the 
ends of the rod by external forces, which we assume are included in an external interaction 
potential U^. 

The angular momentum field has three contributions; from T, U^, and U^, 

L = d,F^ + (l^ + '^0^^ + 2r,^j) + 2g. T'^^- (45) 

i,j,fc=i V i / 1=1 

The functional derivative of if' was evaluated following Eq. but includes an additional 
term derived from the rotation of the frame by variations in q^. There are similar contribu- 
tions from rotations of the frame in the functional derivatives of T and tf . Derivatives of 
the basis vectors di^, and tia with respect to were evaluated using Eqs. IT1.411TT75] from 
Table [B Although the equations of motion must be derived for the canonical momenta Pq, 
and /a, the numerical implementation can use any frame. We have found that it is most 
convenient to use space-fixed linear momenta and body-fixed angular momenta as the pri- 
mary variables, since this seems to minimize the number of rotations of I. The quaternion 
momenta can be rewritten as body-fixed momenta, /j = [tiala + ^iala)/'^-, 

+ E ^^^-^^ = + E '^^^ i^^F^ + r.^f) , (46) 

j,k=l j,k=l 

again using Eq. IT1.5I to evaluate variations in Cia- This expression is equivalent to Eq. [12] 
except that it is written in the body-fixed frame instead of the space-fixed frame. 



B. Discretized Hamiltonian 



In this section we will derive equations of motion for the nodal coordinates and momenta 
by discretizing the line integrals in Eqs. [37] and [iQl The kinetic energy is approximated by 
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the midpoint rule, 



■N 




Yl n 1 ^ iTl 771 \ 



i=l 



where is the discrete kinetic energy per unit length. The discrete Hamiltonian of a set 
of infinitesimal segments, Ti.^ , is an energy density, whereas a Hamiltonian describing finite- 
length segments^i^ would have units of energy. Equation|17]is a second order approximation 
to the kinetic energy of the continuous filament, T = T^As + 0{As)^. Discrete approx- 
imations to the potential energy involve coordinate differences evaluated at the midpoints 
between pairs of nodes. We therefore approximate the potential energy by a trapezoidal 
rule, which is also second order in As, 

^ N 3 

= 2 5Z E^" - - r?) + c?{2e-j: - n^){2eU'' - m ■ (48) 

n=0 i=l 

The derivatives r'^ and q'^ are defined in Eqs. [241 - I25] and the average quaternions q^, used 
to calculate e"„, are defined in Eq. [261 The weights, Wn, for the trapezoidal integration rule 
are w„ = l/2ifn = 0orr;, = A^ and Wn = 1 otherwise. 

The equations of motion for the nodal coordinates and momenta then follow by differen- 
tiation: 

gj-N n 
•n fa 



dpi Mr' ^^^^ 



dl"} 2^ MP 



~ i=i 

/N 



Pi = = (51) 



a 

an-zN ^ f,n inin 

in _ (yrt _ ST^ tj^lj tfc 

i,i,fc=l ^ 

where the nodal forces and torques are 



t 



As 
As 

3 / Q„pC,ra Qn-lpr2,n-l' 



i,j,k=l 



2g" 2q 



r,n-l 



3 / -pn^rjn -pn-l^rjn-l " 



+ E I + ^n-l^r^ ^ ) , (54) 



i,j,k=l 
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and is the length of the unnormahzed quaternion = {q^ + Q'a~^l/2- It is essential that 
the differentiation is done exactly, otherwise the Hamiltonian structure of the equations of 
motion is lost. Equations HHHSSl are straightforward, but Eq. EH requires some explanation. 
The factor of two between the T x and d x contributions (c./. Eq. H^D arises 
because the rate of rotation of the quaternion basis is one-half that of the body-fixed frame. 
Terms involving dot products of with e"^^ vanish by orthogonality, even for the midpoint 
quaternions. Less obviously, the orthogonality of Qa and q'^ is preserved by the discretization, 
so that 

^ (^) ^ 0. (55) 

Although the discrete Hamiltonian, Ti.^ = + U'^ , is only a second-order approximation 
to Ti, the equations of motion for the nodes (Eqs. H9] - [5^ exactly preserve a Hamiltonian 
structure for any As. Equations [331435] do not have this property, although they are the 
same to second order in As. 

For our numerical implementation, it is more convenient to calculate the angular momenta 
in the body-fixed frame rather than the quaternion basis. Making the same transformation 
as from Eq. HS] to Eq. H6l 

3 1 

j,k=l 3 

where the conservative torque in the quaternion basis is given by Eq. [5H No further simpli- 
fication is possible in this case, because the quaternion basis vector is not the same as 
those in the expression for t". The slight variations in the quaternions make the difference 
between the Hamiltonian formulation for the torque (Eq. and the torque (Eq. 1551) derived 
from the finite-difference discretization described in Sec. IIIII 

C. Operator splitting 

Implicit integration methods are typically used to integrate the equations of motion of 
elastic rods,— i^i^ even when the model has no explicit constraints.— The most common 
choice is the implicit midpoint method, which updates the vector Y = [P, Q] to second 
order in the time step At, 

Y{t + At) = Yit) + ^ {Y[Y{t)] + Y[Y{t + At)]) . (57) 
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Implicit methods have the advantage of stabihty for large time steps and the implicit mid- 
point method is in addition symplectic— However a number of force evaluations are needed 
at each time step to solve the non-linear equations (1571) to machine precision, which is nec- 
essary to maintain the symplectic structure. Moreover, the normalization constraint on the 
quaternion is not conserved, 

1=1 

and must be rescaled at each time step. 

Operator splitting techniques are increasingly being used to solve both deterministio^'^'^ 
and stochastic differential equations.—"^ Typically the splitting is devised so that the in- 
dividual propagators can be determined exactly. If the underlying dynamics is strictly 
Hamiltonian,— 1^1^ then symplectic integrators can be constructed by such techniques. The 
Liouville operator, C = C'^ + C'^ , is decomposed into kinetic {C'^) and potential (£^) terms, 

°9r" ' ^'dq' 



ra=l 

here we use a second-order Trotter decomposition,—"^ 

exp [CAt] = exp [C'^At/2] exp [C^ At] exp [C'^At/2] + C(A^)^ (61) 

although higher-order algorithms are available.— 

The integration of the position and momentum equations is a straightforward and exact 
streaming, 

r„(At) = exp [iC^At] r^ = r^ + j^At, (62) 
pl{At) = exp [£^At] = K + f^At, (63) 
/f(At) = exp [C^At] 11 = li + t^At. (64) 

An exact solution of the quaternion update is more complicated, but can be carried out 
using elliptic integrals.— Nevertheless, here we adopt a simpler formulation which uses a 
sequence of rotations about the body-fixed axes, 

n=l \ " 2=1 / I ±a 
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A rotation = l^At/Mp about one of the body-fixed axes changes both the quaternions 
and the other body-fixed momenta: 

exp (£rAt) q: = cos(A0r/2)g„" + sm{A<j>^ /2)C, (66) 

3 

exp (£r At) I] = cos(A0n^ • + E '^^'^ M^<Pm- (67) 

fc=i 

The individual rotations can be combined using any suitable second-order decomposition 
for X^Li "^7^ example 

(exp mAt/2J] exp mAt/2J] exp [£3 At/ J] exp mAt/2J] exp [/:^At/2 J])"^ . (68) 

The update of the quaternions is not exact, but it is symplectic and exactly preserves 
the norm of the quaternion. If the time step is broken up into J subintervals, a more 
accurate integration can be achieved without substantial overhead, since no force evaluation 
is needed.— 

V. NUMERICAL EXAMPLES 

Our analysis has been supplemented by numerical simulations using the algorithms de- 
scribed in the text. We have compared explicit fourth-order Runga-Kutta (RK) integration, 
implicit second-order midpoint (MP) integration, and second-order Operator Splitting (OS) 
(Sec. lIVCp . We have tried each method with forces and torques derived from discretizing 
the partial differential equations (DF), Eqs. ( |Ml) - (l35l) . and with forces and torques derived 
from discretizing the Hamiltonian (DH), Eqs. (I53 i) -( l54l) . We investigated the stability and 
conservation of energy from two initial conditions: a straight filament bent into a circle and 
a straight filament bent into a helix. 

A. A filament bent into a circle 

A straight filament of length 20TTd was bent into a circle of radius lOd and released. The 
dynamics were followed for two different spatial discretizations, dividing the filament into 
63 or 127 equal segments; the corresponding segment lengths were approximately d and 
0.5d. The largest time step for the explicit integrators is Courant limited by the time, tc, 
for a longitudinal wave to cross the shorter of the diameter, d, and the segment length. As; 
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t=300 

t=600 




M 

FIG. 2: Filament shapes at different times: 300to (solid), 600to (long dashes), 900to (dashes), 
1200^0 (dot dash), and ISOOto (dotted). The time scale to = is the time for a longitudinal 
wave to cross the diameter of the filament 

we typically use a time step At = 0.2t(7. As the rod evolves from its initial configuration, 
flexural waves propagate along the filament, leading to a surprising variety of configurations; 
a sampling of the filament shapes is illustrated in Fig. [2l Initially the ends move slowly, 
and the filament assumes a teardrop shape (t = SOOto), followed by a hairpin it = 600to) 
as the ends of the filament accelerate. The time unit to = d/ci, where q is the longitudinal 
wave speed. The inverted U shape (t = 900to) straightens out (t = 1200to), and then 
develops a "double-minimum" shape (t = ISOOto). The center of the filament moves down 
to complete the inversion and the filament approximately retraces the sequence of shapes 
in reverse order, to arrive at the inverted configuration at roughly half the period of the 
main oscillation. However, the motion is not exactly periodic because of the strong coupling 
between the flexural modes. The interaction of flexural waves can lead to large local stresses, 
exceeding that of the initial configuration; for example at the top of the teardrop (t = 300to) 
and at the bends in the hairpin (t = 600to). It has been shown that flexural modes can 
cause unexpected fractures by this mechanism.— 

A complete cycle of the filament motion, back to a rough approximation of its initial 
configuration, takes about 6000to for a filament of length L ~ 60d, and is quadratic in 
the length of the filament. The scaling is due to the dispersion relation of flexural waves, 
uj oc k"^, which is quadratic rather than linear in the wavevector (k); the period of the longest 
flexural wave. Sir/ {cik'^d) is roughly lO^to. A plot of energy vs. time. Fig. [3^, shows that all 
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FIG. 3: Conservation of energy for symplectic (OSDH) and non-symplectic (RKDF, RKDH, OSDF) 
algorithms. The initially circular configuration of the filament unwinds as illustrated by the snap- 
shots in Fig. [21 a) 63 segments, At = 0.2to = 0.2tc; b) 63 segments, At = 0.02to = 0.02tc; 
c) 127 segments, At = O.Uq = 0.2tc', d) OSDH algorithm with varying precision, 63 segments, 
At = 0.2to = 0.2tc. 



the algorithms integrate stably for about 10 oscillations, but only the symplectic methods, 
MPDH and OSDH, are stable at long times; on the scale of Fig. [31 results for MPDH and 
OSDH superpose, so only the results for OSDH are shown. We have run the MPDH and 
OSDH algorithms to a time of lO^to or 16000 periods, with no indication of instability. By 
contrast, changing the forces to the non-Hamiltonian form (OSDF) or switching to the RK4 
integrator (RKDH) causes instabilities at times of the order of lO^to- Reducing the time 
step. Fig. [3)3, improves the stability of the Runga-Kutta integration of the Hamiltonian 
forces (RKDH), increasing the range of stability by about an order of magnitude. This is 
because RKDH becomes symplectic in the limit At 0. On the other hand if the forces 
are not Hamiltonian, reducing the time step does not improve the stability; both RKDF 
and OSDF algorithms become unstable after a time of about lO^to, regardless of time step. 
The discretized forces approach a Hamiltonian form in the limit As — and reducing 
the segment length improves the stability of the OSDF algorithm, extending the range of 
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stability by about a factor of 4 for a twofold reduction in the segment length, Fig. [St. 
However, this is a double limiting process requiring a progressively smaller time step as well 
as a reduced segment length, making it computationally expensive. The RKDF algorithm 
is not helped by a reduction in segment length; it needs a further reduction in time step as 
well to see any improvement. 

The non-linearity of the dynamics causes the filament to eventually reach a state of ther- 
mal equilibrium, fluctuating around the straight configuration. For the 63 segment rod the 
equilibration time is about lO'^to independent of time step. For a constant filament length, 
we observe that the equilibration time is roughly quadratic in the number of segments. Thus 
the behavior of this system in the continuum limit is an interesting question for future work, 
but beyond the scope of the present paper. 

The stability of the symplectic integrator is affected by accumulated round-off error. The 
results in Fig. [3li show that the symplectic integration scheme (OSDH) is quite unstable in 
single precision arithmetic. The most rapid instability, at t < lO'^to, "was traced to accu- 
mulated errors in the quaternion normalization. The operator splitting algorithm maintains 
the quaternion normalization to machine precision and with 64-bit arithmetic the normaliza- 
tion error is stable at less than one part in 10^^. But in single precision, the error increases 
rapidly, which causes an incompatibility with the assumption that the nodal quaternions are 
normalized. More puzzling is that rescaling the quaternions does not solve the problem, but 
merely delays the onset of the instability. However, if the initial accumulation of round-off 
error is random, we would expect the double precision version to run stably for about 10^^ 
times longer, or lO^^to which is well beyond the event horizon of the simulation. 

The short-time fluctuations in energy of the OSDH algorithm cannot be seen on the scale 
of Fig. [3l but they are quadratic in the time step, with a relative magnitude of approximately 
0.1(At/to)^- These short-time fluctuations in energy are about 20 times larger with OSDH 
than with MPDH. However there is also a drift in the energy with time, again quadratic 
in At, but larger, as shown in Fig. |H Over long time intervals, OSDH preserves energy 
conservation with about an order of magnitude better accuracy than MPDH at the same At 
(Fig. Hj). MPDH requires 5-10 times as many force evaluations as OSDH per time step, so 
that the explicit operator splitting algorithm is clearly preferable for long-time dynamics. 

Dichmann and Maddocks studied the dynamics of a Kirchoff rod from the same initial 
configuration,— but with the filament pinned at one end. The nodal forces and torques 
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FIG. 4: Conservation of energy for symplectic algorithms OSDH and MPDH; 63 segments were 
used in each case, a) OSDH, At = 0.2to; b) MPDH At = 0.2to; c) OSDH, At = 0.02to; d) MPDH, 
At = 0.02to. 

were also Hamiltonian, but the implicit midpoint integrator was used instead of operator 
splitting. Their results showed a small drift in the total energy of around 0.2% after approx- 
imately 30 oscillations of the filament, or 200, OOOto in our units. Our results for the MPDH 
algorithm behave in a qualitatively similar fashion; with a time step At = 0.2to we observe 
an accumulated energy drift of 0.3% at t = 200, OOOto- The error with OSDH is about an 
order of magnitude smaller. The GE model requires a smaller time step to explicitly inte- 
grate the shear and extensional degrees of freedom, but surprisingly, it is only a factor of 
8 smaller than the time step used for the constrained rod.— This suggests that the explicit 
OSDH algorithm can integrate the full GE rod model with about the same computational 
cost as an implicit integration of the Kirchoff model. If excluded volume interactions are 
included, it is likely that these very stiff forces will set the overall time step, as is typical 
in molecular dynamics simulations. In such cases the computational advantages of a fully 
explicit simulation will be considerable. 
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a) b) 




FIG. 5: Filament shapes at different times: a) t = 0; b) t = lOOto; c) t = 200to; d) t = SOOto; e) 
t = 400to; f) i = SOOto- The simulations with 630 segments are shown as thick solid lines, while 
simulations with 63 segments are shown by the spheres. 

B. A filament bent into a helix 

We have also examined a more complicated initial condition, a straight rod of length 
207r(i wound into a tight helix with exactly four complete turns. The curvature, f2 = 
[0.4(i~^, 0, O.lc?"^], is high and generates motion in all three spatial dimensions, which poses 
a difficult challenge for the numerical method. We used two different discretizations, 63 
segments of length As ~ d and 630 segments of length As ^ O.ld; snapshots of the initial 
evolution of the filament shapes are shown in Fig. O There is a high degree of dynamical 
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FIG. 6: Conservation of energy and thermal equilibrium with the symplectic integrator OSDH. 
The initially helical configuration of the filament unwinds as illustrated by the snapshots in Fig. [3 
The kinetic, potential and total energy of the 63 segment model (At = O.lto = O.ltc) are shown 
for: a) lO^o and b) 4000to. The body-fixed kinetic energy of the individual degrees of freedom is 
also shown: c) Shear and extension and d) Bending and torsion. 

coherence between the results at the two different resolutions, although the strong nonlin- 
earity of the problem means that they start to diverge at times of the order of SOOto- We 
did not include any excluded volume interactions in these simulations, and the filaments can 
therefore cross; this does not affect the accuracy of the numerical algorithm. 

As in the planar bend case, the symplectic algorithm (OSDH) conserves energy. Fig. [6^, 
for as long a time as we have tested, up to lO^to- The non-linear coupling is much stronger 
than in the previous example, because of the higher curvature and the three-dimensional 
deformation; here the filament rapidly comes to thermal equilibrium. The loss of coherent 
oscillations can be seen more clearly in the expanded time scale of Fig. [6}3. Over the same 
time scale, lO'^to, we see that equipartition of energy is established between the various 
degrees of freedom. Figs. [Ht and [Hli; similar results holds for the various components of the 
potential energy as well. Unlike the planar bend case, here the more finely resolved filament 
(630 segments) comes to thermal equilibrium on more or less the same time scale, ~ 40, OOOto, 
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rather than lO^to as would be expected for a quadratic scahng of the equihbration time with 
A^. This suggests fundamental differences in the dynamics of the two-dimensional bending 
from the full three-dimensional problem. 

VI. CONCLUSIONS 

In this paper we have presented a new algorithm for simulating the dynamics of elastic 
filaments. The test problems show the method to be extremely stable, with exact con- 
servation of momentum and angular momentum (to machine precision), and global energy 
conservation to order At^. The algorithm is fully explicit and requires no constraints of any 
kind, neither on the forces nor on the quaternions. It is thus simpler in some ways than 
typical WLC implementations which include extensional forces as a constraint. In contrast 
to the WLC, the GE model correctly incorporates large bending deformations and twisting; 
it includes the Kirchoff rod as a limiting case. 

Symplectic integration of the GE model can use a large time step, within a factor of 10 
of a constrained filament^ that excludes shear and extensional modes. Explicit operator 
splitting has better long-term energy conservation than the implicit midpoint method and 
requires an order of magnitude fewer force evaluations per time step. In cases where the 
time step is limited by the stiffness of excluded volume interactions, the GE model may be 
more computationally efficient than the Kirchoff model, due to the absence of constraints. 

In this work we only discussed Hamiltonian systems, but operator splitting is a powerful 
method for integrating stochastic systems as well.— We have considered the case when the 
rod is subjected to dissipative and random forces, in addition to the elastic forces. Using 
operator splitting we can integrate the momentum equation exactly, using the Ornstein- 
Uhlenbeck solution, and therefore preserve quadratic norms to order At^, as opposed to the 
At accuracy of Brownian dynamics. This work will be reported in a future paper. 
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APPENDIX A: PROPERTIES OF QUATERNIONS 



A quaternion Z = qQ + qxi + Qyj + qzk is a complex number with multiplicative identities 

i"^ = f = = ijk = -1. (Al) 

We use the notation Z to indicate the quaternion and to denote a vector containing the 
scalar, go, and vector, q = [qx,(ly,(lz], components of Z. The quaternion algebra, Eq. flAll) . 
leads to rules for multiplication that are analogous to the cross-product of unit vectors: 

ij = -ji = k 

jk = -kj = I (A2) 
ki = —ik = j 

If we then identify i, j, k, with Cartesian unit vectors i, j, k, the multiplication of two 
quaternions, Z = go + q^i + qyj + qzk and Z' = qQ + q'^i + q'yj + g^fc, can be written, using 
the quaternion multiplication rules defined in Eqs. lAll and lA2t as 

ZQZ' = qoq'o - q ■ q' + qoq' + q'oq + q x q\ (A3) 

where denotes a quaternion multiplication. 

A vector u can be rotated by the unitary transformation Z QuQ Z^^, where the multi- 
plicative inverse of a unit quaternion is Z~^ = Qo — Q- Applying Eq. IA3I and treating u as 
a quaternion with zero scalar component, the rotated vector u' is given by 

u' = {ql- q - q)u + 2qq ■ u + 2gog x it, (A4) 

and remains a pure vector. The rotation can also be written in matrix form, u[ = diaUa, 
with the director vectors that form the rotation matrix dia as given in Eq. IT1.2I 
An infinitesimal change in the directors is given by a rotation 6(p: 

6di = 6(f) X di, (A5) 

with ft = ds4> and u: = dt4>. The combination of the original rotation Z and an additional 
infinitesimal rotation 6Z = 1 -|- 6(p/2 can be found by applying the rotations sequentially, 

u' + 6u = 6Z Q Z QuQ Z~^ Q 6Z-^ = Z' Q u Q Z'-\ (A6) 
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The quaternion Z' is found by multiplying the two quaternions, 



Z' = bZ(dZ = qQ-q--^^q^ qo— q x 



6(f) 
T 



6(f) 

Y 



6(t) 



(AT) 



Thus, the variation in the quaternion 6Z = Z' — Z is linearly related to 6(f>, 

f ^<lo\ ( —o-r —a., —a. \ 
6qx 
6qy 



(A8) 



-qz Qo Qx 
K^Qz J \ Qy -Qx qo ) 

The column vectors in Eq. IA8I define a set of basis vectors in the quaternion space, e^a, 
where taa is the transpose of the matrix in Eq. IA8[ These basis vectors are orthogonal to 
qa and relate changes in quaternions to rotations about the space-fixed axes. 



2e„a^'?a, 6qa = -eaa6(f)a- 



(A9) 

In this work we have used body-fixed rotations, Eqs. [2T| - [22l for which we need the basis 
vectors e^a given in Eq. IT1.3[ they are related to the space fixed basis Caa by the rotation 
matrix, e^a = diaCaa- The vectors 6^ or e^a, together with g^, form a complete basis in the 
quaternion space. 

Finally, we obtain the derivatives of the basis vectors quoted in Eqs. IT1.4HTL51 A 
variation in the basis vectors di is related to an infinitesimal rotation, Eq. IA5t 

3 3 

6dia = €ai3-y6(p/sdi^ = ^ijkdja64>k = 2 ^ijkdja^kb^qb- (^10) 

j,k=l j,k=l 

The variation in di can also be directly related to constrained variations in quaternions, 

dd^ 



6di, 



dqa 



i6ab - qaQb) 6qb 



(All) 



where the projection operator {6ab — QaQb) is included to ensure that the normalization con- 
dition, 6qaqa = 0, is satisfied. Equation IT1.4I can then be obtained by making use of the 
result 



ddi, 



dqa 



2d- 



(A12) 



The rotation matrix can be written as a product of e vectors, di^ = eiaCaa- A space 
fixed vector is first rotated into the quaternion basis by Caaf^ and then rotated from the 
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quaternion basis to the body- fixed frame by 2eia- A variation in dia is tlien composed of two 
equal contributions from variations in e^a and Cia, 

Substituting Eq. flAlOp for the variation in dia-, and using the orthogonahty of the d vectors, 

3 
k=l 

Multiplying both sides by Cji, and summing over j, 

3 

{^ab - Qaqb) Seia = ^ (^ijkejbeka^qa- (A15) 

j,k=l 

The variation in can also be related to constrained variations in g^, c./. Eq. (lAlip . using 
the relation qaScia = -eia6qa, 

3 

- (iSi.. — nunS\ iSrii. = 

dqc 

Equation IT1.5I then follows from 



Seia = {6bc - qbqc) Sqb = 2J (^ijkCjaekbSqb - qa^ib^qb- (A16) 

j,k=i 



Oc ■ 

qb^^ = Gia- (A17) 

oqb 
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